home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / FIXRTS.DEM < prev    next >
Text File  |  1991-04-29  |  2KB  |  81 lines
  1. PROGRAM d12r11(input,output);
  2. (* driver for routine FIXRTS *)
  3. CONST
  4.    npoles=6;
  5.    twonp=12;   (* twonp=2*npoles *)
  6.    twonp2=14;   (* twonp2=twonp+2 *)
  7. TYPE
  8.    gl2array = ARRAY [1..2] OF real;
  9.    glnparray = ARRAY [1..npoles] OF real;
  10.    glcarray = ARRAY [1..twonp2] OF real;
  11. VAR
  12.    i,j : integer;
  13.    dum : real;
  14.    polish : boolean;
  15.    d : glnparray;
  16.    sixth : gl2array;
  17.    zcoef,zeros : glcarray;
  18.  
  19. (*$I MODFILE.PAS *)
  20. (*$I LAGUER.PAS *)
  21.  
  22. (*$I ZROOTS.PAS *)
  23.  
  24. (*$I FIXRTS.PAS *)
  25.  
  26. BEGIN
  27.    d[1] := 6.0; d[2] := -15.0; d[3] := 20.0;
  28.    d[4] := -15.0; d[5] := 6.0; d[6] := 0.0;
  29.    polish := true;
  30.    (* finding roots of (z-1.0)^6 := 1.0 *)
  31.    (* first write roots *)
  32.    zcoef[2*npoles+1] := 1.0;
  33.    zcoef[2*npoles+2] := 0.0;
  34.    FOR i := npoles DOWNTO 1 DO BEGIN
  35.       zcoef[2*i-1] := -d[npoles+1-i];
  36.       zcoef[2*i] := 0.0
  37.    END;
  38.    zroots(zcoef,npoles,zeros,polish);
  39.    writeln('Roots of (z-1.0)^6 = 1.0');
  40.    writeln('Root':22,'(z-1.0)^6':27);
  41.    FOR i := 1 to npoles DO BEGIN
  42.       sixth[1] := 1.0;
  43.       sixth[2] := 0.0;
  44.       FOR j := 1 to 6 DO BEGIN
  45.          dum := sixth[1];
  46.          sixth[1] := sixth[1]*(zeros[2*i-1]-1.0)
  47.             -sixth[2]*zeros[2*i];
  48.          sixth[2] := dum*zeros[2*i]
  49.             +sixth[2]*(zeros[2*i-1]-1.0)
  50.       END;
  51.       writeln(i:6,zeros[2*i-1]:12:6,zeros[2*i]:12:6,
  52.          sixth[1]:12:6,sixth[2]:12:6)
  53.    END;
  54.    (* now fix them to lie within unit circle *)
  55.    fixrts(d,npoles);
  56.    (* check results *)
  57.    zcoef[2*npoles+1] := 1.0;
  58.    zcoef[2*npoles+2] := 0.0;
  59.    FOR i := npoles DOWNTO 1 DO BEGIN
  60.       zcoef[2*i-1] := -d[npoles+1-i];
  61.       zcoef[2*i] := 0.0
  62.    END;
  63.    zroots(zcoef,npoles,zeros,polish);
  64.    writeln;
  65.    writeln('Roots reflected in unit circle');
  66.    writeln('Root':22,'(z-1.0)^6':27);
  67.    FOR i := 1 to npoles DO BEGIN
  68.       sixth[1] := 1.0;
  69.       sixth[2] := 0.0;
  70.       FOR j := 1 to 6 DO BEGIN
  71.          dum := sixth[1];
  72.          sixth[1] := sixth[1]*(zeros[2*i-1]-1.0)
  73.             -sixth[2]*zeros[2*i];
  74.          sixth[2] := dum*zeros[2*i]
  75.             +sixth[2]*(zeros[2*i-1]-1.0)
  76.       END;
  77.       writeln(i:6,zeros[2*i-1]:12:6,zeros[2*i]:12:6,
  78.          sixth[1]:12:6,sixth[2]:12:6)
  79.    END
  80. END.
  81.